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Abstract. We present a 3D N-body simulation of the Milky Way in- 
cluding 4-10^ star- and dark-like particles and 2-10'^ gas particles, initially 
distributed according to an axisymmetric, observationally constrained, 
mass model. The whole system is self-gravitating and the gas hydrody- 
namics is solved using the SPH method. The simulation leads to the 
spontaneous formation of a central bar that strongly affects the gas dy- 
namics. We compute [l — V) diagrams for both the gaseous and the stellar 
particles as a function of the angle of the bar with respect to the observer 
and compare the results for the gas with HI and CO observations. 



1. Introduction 

Since the first direct evidences for the Galactic bar (Blitz & Spergel 1991; Wein- 
berg 1992), gas and stellar dynamics in the Milky Way have experienced a clear 
renewal of interest. Several new photometric and kinematic data have already 
been or are on the way to be used to constrain analytical triaxial bulge models 
(e.g. Stanek, this volume). By performing numerical simulations, the prob- 
lem can be tackled differently: time sampled snapshots of the numerical galaxy 
can be compared to the observations varying various model parameters (see e.g. 
Wada et al. 1994). Invaluable informations can then be inferred concerning 
quantities hard to observe, like the bar pattern speed. Nevertheless, so far none 
of the dynamical models have included in a completely self-consistent way three 
dimensional gaseous, stellar, and dark components. This is however a prerequi- 
site to highlight various global asymmetries or instabilities, or if the evolution 
of the model needs to be followed for Gyrs since secular evolution processes are 
especially active in barred galaxies (see e.g. Martinet 1995). 



2. Initial Conditions 

Our initial axisymmetric 3D mass distribution includes four components (see Ta- 
ble 1): an exponential stellar disk of constant thickness, a composite power-law 
stellar nucleus-spheroid (NS), a dark halo (DH) to ensure an almost flat rotation 
curve out to about 30 kpc, and a dissipative gas component consisting of a kind 
of smoothly truncated Mestel disk with a central core and an outward linearly 
increasing scale height. The gas mass inside the solar circle is 4.1 • 10^ Mq. The 
parameters of this mass model have been adjusted to several observational con- 
straints, like the rotation curve and local mass densities, assuming Ro = 8 kpc. 
To each of the stellar and dark components is assigned an individual isotropic 



1 



Table 1. Initial mass model of the simulation. The oblate NS and DH 
have a common axis ratio e = 0.5. 



Components # particles Space density Parameters 

Disk 150, 000 cx exp {-R/hn) ■ sech {z/h^,) Hr = 2.5 kpc 

= 250 pc 



^0 



NS 50,000 cx , P = -1.8, g = -3.3 

1 + "^^"^ a = 1 kpc 

= + z^le'^)/a^ M = 3.0 • 10^° 

DH 200,000 (xexp(-/i), 6 = 9.4 kpc 

H'^ = {B? + z''le^)l\? M = 3.1 • IQii 



Gas 20,000 cx exp (-||^ - i^l^) s = 0.017 



Rg = 20 kpc 



M = 1.13 • 10^° Mr 



Maxwellian velocity distribution such that the velocity dispersion and the mean 
azimuthal velocity (no streaming motion in other directions) satisfy the Jeans 
equations. Initially, the gas particles have circular velocities and a uniform 
sound speed of 10 km/s, providing only a rough hydrostatic equilibrium. The 
non-inclusion of star formation presently is the strongest limitation of the model. 



3. Time Evolution 

The system has been integrated for 2.6 Gyr using the PM method on polar- 
cylindrical grid to compute the gravitational forces (Pfenniger &; Friedli 1993) 
and the Lagrangian SPH technique (see e.g. Benz 1990) with isothermal internal 
energy to compute the pressure and viscous forces acting on the gas. The central 
radial resolution of the grid is 38 pc and the gas smoothing length near the center 
is of the same order. The time integrator is a synchronized leap-frog with an 
adaptativc time-step At < 0.1 Myr able to fully resolve shocks (variable At in 
time but constant for all particles). The evolution is completely self- consistent 
and free from any imposed symmetry. 

The simulation, illustrated in Fig. 1, leads to a growing central bar after 
about 1 Gyr, while m = l modes appear in the stellar disk. The bar itself never 
really behaves like a bi-symmetric solid rotating body, but instead undergoes 
many asymmetrical distortions, particularly strong during its first rotations. 
Such perturbations also induce irregular gas flows (like the one-side void inside 
3 kpc occurring at t = 1.5 Gyr) capable of producing peculiar signatures in 
the gaseous {I — V) diagrams. The gas particles near the center condense into 
a nuclear ring which progressively shrinks and finally forms a small disk of 
about 300 pc radius. In this process, the ring suffers several instabilities, like 
tilting, warping and off-centering. There is also an increasing central gas core 
accreting gas from the nuclear ring/disk. The final properties of the bar are a 
semi-major axis a « 3.5 kpc, an axis ratio h/a ^ 0.6 — 0.7 and a pattern speed 
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Figure 1. Face-on configuration of the visible components at various times. 
Left: stellar (disk+NS) isodensity curves. Middle: gas particles, same scale. 
The position of the observer in Figs. 2d and 3b is indicated by a solid circle. 
Right: zoom on the central gas. The encircled dots at 1.5 Gyr show the par- 
ticles responsible for the elongated strip in Fig. 3b. The length unit is 1 kpc. 
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rjp Sri 50 km/s, which places the CR around 4 kpc, the ILR shghtly outside 1 kpc 
and the OLR close to 7 kpc. A comparison with symmetrized simulations (i.e. 
imposed reflexion symmetries about the z-axis and the z = plane) without gas 
(Fux et al. 1995) and with gas indicates that the imposed symmetries generate 
more extended bars, whereas the presence of gas makes them slightly rounder. 
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Figure 2. a) Definition of the angle (/? of the bar relative to the ob- 
server, if is positive in the direction of rotation, b) HI (/ — V) map for 
< / < 220° (Hartman & Burton 1995). c) CO {I - V) map (Dame 
et al. 1987). d) Distance weighted (l — V) diagram from the simulation 
at i = 1.7 Gyr and for ip ^ 30°. All (l — V) plots are averaged over 
|6|<2.25°. Higher density regions are whiter. 



4. Gaseous (/ — V) Diagrams 

For the moment, the longitude- velocity diagrams of the gas particles have only 
been explored every 100 Myr, as a function of the relative angle if of the bar with 
respect to the observer located at R = Ro (see Fig. 2a). It comes out that for 
a given angle cp, these diagrams are time-dependent. Moreover, for each value 
of if there is a modulo tt choice in orienting the bar, allowing for two possible 
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Figure 3. a) Details of the CO (l — V) map in the Galactic central 
region. Note the positive high velocity emission strip reaching very 
negative longitudes, b) One of the most resembling gaseous {I — V) 
plot available from the simulation at t = 1.5 Gyr and with (p ~ 20°, 
restricted to the gas particles inside |6| < 2.25°. 



[l — V) diagrams. These would look identical under reflexion symmetry about 
the z-axis, but can appear very different in our symmetry-free case. 

When computing (I — V) diagrams over a large longitude range in order 
to confront them with observations, one should weight the contribution of each 
particle by its inverse squared distance relative to the observer. Figure 2 shows 
such a diagram for a realistic position angle of the bar at a time where similarities 
with the observed HI and CO maps clearly emerge. The gas in our simulation 
can be considered as a mixture of atomic and molecular gas. 

A well-known puzzling feature in the observed {I — V) data around the 
Galactic center is the presence of gas in regions forbidden to pure circular motion. 
Part of it can be explained by the signature of bar elongated xi or anti-bar X2 
periodic orbits (Binney et al. 1991) or by more sophisticated gas flow modeling 
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(Mulder & Liem 1986; Jenkins & Binncy 1994) in a fixed bar potential. However, 
these approaches do not account for the intriguing gas strip at F < 100 km/s 
extending in the negative longitude region down to about I = —5° (see Fig. 3a). 
Our simulation displays at least once such a peculiarity, as shown in Fig. 3b. 
The spatial positions of the gas particles populating the "forbidden" strip in 
this figure are highlighted in Fig. 1. A closer look to the simulation reveals 
that these particles were previously strongly shocked near the bar end, then 
quickly fell toward the center to finally overtake it with high positive radial 
velocities. Clearly, the high degree of asymmetry displayed by the observations 
and achieved here, for instance at t = 1.5 Gyr, cannot occur in standard bi- 
symmetric bar models. 

5. Stellar {I — V) Diagrams 

Constraining the bar properties using stellar (I — V) diagrams will probably 
prove to be at least as powerful as the methods based on gas dynamics. Since 
observations in this field mainly rely on star counts, corresponding model dia- 
grams do not have to take into account a distance weighting factor as for the gas. 
Fig. 4 gives a first example of what can be expected from such diagrams. The 
signatures obviously depend on the bar orientation. A detailed confrontation 
with observations will be published elsewhere. 
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Figure 4. Distance unweighted {I — V) diagram of the stellar disk as 
a function of (p at t = 2.5 Gyr, averaged over all b. Diagrams for (p and 
if + Tr have been added together to reduce noise. 
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6. Conclusion 



A fully consistent numerical 3D barred model of the Milky Way in good qualita- 
tive agreement with HI and CO kinematic observations at both small and large 
longitudes has been presented. In particular, morphological and kinematical 
asymmetries naturally occur during the simulation which are able to reproduce 
the high velocity CO feature at negative longitudes in the (Z — V) diagram. 

Discussion 

Christodoulou: Your simulation run for just a few Gyrs. If you managed to take 
it to 15 Gyr, which of the features would you expect to survive (bar, spirals, 
nuclear ring)? How much mass is accumulated in the nuclear ring? 

Fux: Integration over a Hubble time would be unrealistic because the initial 
conditions mimic the present state of the Milky Way and the absence of star 

formation favor the growing of the central gas concentration, thus shortening 
the life-time of the bar. At 1.5 Gyr, the ring mass is about 6 • 10^ M0. 

Stanek: Is your evolved disk similar to the real disk of the Galaxy (i.e. scale 
height, velocity ellipsoid, etc.)? 

Fux: At t = 1.5 Gyr, the scale height of the stellar disk is about 350 pc outside 
corotation and its maximum l.o.s. velocity dispersion towards Baade's Window 
about 130 km/s. Its velocity ellipsoid has become spontaneously anisotropic with 
local velocity dispersions close to the observed ones for the Galactic old disk. 

Teuben: To what extent are the features like a tilted or offset nuclear disk due 
to \/]V-noise, or robust. Is the tilt related to vertical instability strip? 

Fux: They are robust with a density contrast far above the one resulting from 
noise fluctuations. A vertical instability strip does exist, but we haven't yet 
checked that the associated orbit families are consistent with the tilt morphology. 
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